Scaling of stiffness energy for 3d ±J Ising spin glasses 
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Large numbers of ground states of 3d EA Ising spin glasses 
are calculated for sizes up to 10"^ using a combination of a 
genetic algorithm and Cluster-Exact Approximation. A de- 
tailed analysis shows that true ground states are obtained. 
The ground state stiffness (or domain wall) energy A is cal- 
culated. A |A| ~ L®s behavior with Qs = 0.19(2) is found 
which strongly indicates that the 3d model has an equilibrium 
spin-glass-paramagnet transition for non-zero Tc- 

Keywords (PACS-codes): Spin glasses and other 
random models (75.10.Nr), Numerical simulation studies 
(75.40.Mg), General mathematical systems (02.10.Jf). 

Introduction The question whether three dimensional 
(3d) Ising spin glasses have a non-zero transition tem- 
perature Tc is still not answered beyond all doubts. The 
best evidence for ordering below a finite Tc was found 
recently [^,|| by extensive Monte Carlo Simulations. But 
the authors could not completely rule out other scenarios. 

In this work we address the question by calculating 
the stiffness or domain wall energy A = — E'p which 
is the difference between the ground state energies E"-, 
for antiperiodic and periodic boundary conditions 
This quantity was studied earlier only for very small sys- 
tem sizes using T = transfer matrix methods and 
Monte Carlo simulations The stiffness energy shows 
a finite size dependence 
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where L is the linear system size. A positive value of the 
stiffness exponent indicates the existence of a spin 
glass phase for non-zero temperature. Since the direct 
calculation of ground states for 3d spin glasses is NP- 
hard there is no polynomial algorithm available. In our 
work we calculate ground states using a combination of 
Cluster-Exact Approximation (CEA) and a genetic 
algorithm Similar calculations proved that the 2d 

spin glass, where exact ground states can be calculated 
using a polynomial time algorithm, exhibits no ordering 
for T > 0. 

We investigate systems of N spins ai — ±1, described 
by the Hamiltonian 
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neighbor interactions (bonds) take independently — 
±1 with equal probability. The antiperiodic boundary 
conditions for calculating are realized by inverting 
one plane of bonds. 

Algorithm The algorithm for the calculation bases on 
a special genetic algorithm and on Cluster-Exact 
Approximation ^ which is an sophisticated optimiza- 
tion method. Now a short sketch of these algorithms is 
given, because later the influence of different simulation 
parameters on the results is discussed. 

The genetic algorithm starts with an initial popula- 
tion of Mi randomly initialized spin configurations {— 
individuals), which are linearly arranged in a ring. Then 
i/ X Mi times two neighbors from the population are taken 
(called parents) and two offsprings are created using a tri- 
adic crossover: a mask is used which is a third randomly 
chosen (usually distant) member of the population with a 
fraction of 0.1 of its spins reversed. In a first step the off- 
springs are created as copies of the parents. Then those 
spins are selected, where the orientations of the first par- 
ent and the mask agree The values of these spins 
are swapped between the two offsprings. Then a muta- 
tion with a rate of pm is applied to each offspring, i.e. a 
fraction pm of the spins is reversed. 

Next for both offsprings the energy is reduced by ap- 
plying CEA. The method constructs iteratively and ran- 
domly a non-frustrated cluster of spins, whereas spins 
with many unsatisfied bonds are more likely to be added 
to the cluster. For 3d ± J spin glasses each cluster con- 
tains typically 58 percent of all spins. The non-cluster 
spins act like local magnetic fields on the cluster spins. 
For the spins of the cluster an energetic minimum state 
can be calculated in polynomial time by using graph the- 
oretical methods |12|-14|: an equivalent network is con- 
structed the maximum flow is calculated |l^ , [l7| p] 
and the spins of the cluster are set to their orientations 
leading to a minimum in energy. This minimization step 



In this letter we consider 3d cubic systems with peri- 
odic boundary conditions, N = spins and the nearest 



^Implementation details: We used Tarjan's wave algorithm 
together with the heuristic speed-ups of Traff. In the con- 
struction of the level graph we allowed not only edges (w, w) 
with level(u') = level (u) -1-1, but also all edges {v,t) where t is 
the sink. For this measure, we observed an additional speed- 
up of roughly factor 2 for the systems we calculated. 
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is performed n^nin times for each offspring. 

Afterwards each offspring is compared with one of its 
parents. The pairs are chosen in the way that the sum 
of the phenotypic differences between them is minimal. 
The phenotypic difference is defined here as the num- 
ber of spin positions where the two configurations differ. 
Each parent is replaced if its energy is not lower (i.e. bet- 
ter) than the corresponding offspring. After this whole 
step is done v x Mi times the population is halved: From 
each pair of neighbors the configuration which has the 
higher energy is eliminated. If not more than 4 individu- 
als remain the process is stopped and the best individual 
is taken as result of the calculation. 

The following representation summarizes the algo- 
rithm. 

algorithm genetic CEA({Jy }, Mi, p„, nmin) 
begin 

create Mi configurations randomly 

while [Mi > 4) do 

begin 

for i = 1 to ly X Mi do 

begin 

select two neighbors 

create two offsprings using triadic crossover 
do mutations with rate pm 
for both offsprings do 
begin 

for j = 1 to riniin do 

begin 

construct unfrustratcd cluster of spins 
construct equivalent network 
calculate maximum flow 
construct minimum cut 
set orientations of cluster spins 
end 

if offspring is not worse than related parent 
then 

replace parent with offspring 

end 
end 

half population; Mi — Mi/ 2 
end 

return one configuration with lowest energy 
end 



The whole algorithm is performed nji times and all 
configurations which exhibit the lowest energy are stored, 
resulting in Hg statistical independent ground state con- 
figurations. 

This algorithm was already applied to examine the 
ground state structure of 3d spin glasses Q . 

Results For each system size we tried many dif- 
ferent combinations of the simulation parameters 



mi,v,nmin,Pm for some sample systems. The final pa- 
rameters where determined in a way, that by using four 
times the numerical effort no reduction in energy was 
obtained. Here Pm = 0.2 and tir = 10 were used for all 
system sizes. Table | summarizes the parameters. Also 
the typical computer time r per ground state computa- 
tion on a 80 MHz PPC601 is given. 

Ground states were calculated for system sizes up to 
i = 10 for Nl independent realizations (see table |) of 
the random variables. For each realization the ground 
states with periodic and antiperiodic boundary condition 
were calculated. One can extract from the table that the 
L = 10 systems alone required 1990 CPU-days. Using 
these parameters on average Ug > 8 ground states were 
obtained for every system size L using ur = 10 runs per 
realization. 

In the second part of this paragraph a detailed analysis 
of the influence of the simulation parameters is presented. 
But at first the results for the stiffness energy are shown 
in Fig. g. 

Also given is a fit A(i) - which results in 65 = 
0.19(2). Because of the large sample sizes the error bars 
are small enough, so we can be pretty sure that Qs > 0. 
It means that the 3d EA spin glass exhibits a non-zero 
transition temperature Tc- Since A is a small difference 
of large values, we have to be sure that we calculate true 
ground states in order to believe our results. 

Fig. H shows the average energy per spin for 90 test 
systems of size L = 6 as function of the product of the 
simulation parameters Mi x v x n^i^ which is propor- 
tional to the computer time since all other parameters 
are kept fixed. The sets {Mi, i^,n^in) = (8, 1, 1), (8, 2, 1), 
(8,4,1) and M, = 8,16,32,64,128 for (i^, n^in) = (4,2) 
were investigated. The energy decreases with increasing 
numerical effort. For the sets with Ali > 32 the energy 
does not decrease further. We compared our results for 
the test systems with exact ground states which were ob- 
tained using a Branch-and-Cut program |]l9| , ^ . For the 
two largest parameter sets the genetic CEA algorithm 
found the true ground states for all 90 systems! The 
same result was obtained for L = 4 as well.|^ So we can 
be sure that genetic CEA and our method of choosing 
the parameters lead to true ground states or at least to 
states very close to true ground states. 

The choice of the simulation parameters has only a 
small influence on the the energy for large values of M^. 
A more sensitive indicator is the stiffness energy as func- 
tion of the simulation parameters. This is shown for the 
same 90 test systems in Fig. |[ We also calculated the 
exact ground states for the realizations with antiperiodic 
boundary conditions and found again that genetic CEA 



^For L > 6 the Branch-and-Cut needs to much computer 
time because of the exponential time complexity. 
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produced the exact results for the parameter sets with 
M, > 32. 

Since the parameters were tested always on a restricted 
number of systems we can not absolutely be sure that 
genetic CEA always finds true ground states. Since we 
are interested in the stiffness energy A we take a closer 
look at it. If states very close to the true ground states 
are found, the resulting stiffness energy may be for some 
realizations smaller and for others higher than the correct 
result. So we expect that the effect should cancel out with 
increasing sample size. 

This is confirmed by Fig. |^ where the stiffness energy 
as function of the sample size Nl is shown for three dif- 
ferent parameter sets (8,4,1), (16,4,2) and (64,4,2). The 
stiffness energies found using the second set, where we 
are very close to the true ground states, converges to the 
values found using the third set. The values found us- 
ing the first set do not converge. So even by calculating 
states very close to the true ground states on gets very 
good estimates of the stiffness energy. 

Conclusion Results have been presented from calcu- 
lations of a large number of ground states of 3d Ising 
spin glasses. They were obtained using a combination of 
Cluster-Exact Approximation and a genetic algorithm. 
The finite size behavior of the stiffness energy has been 
investigated. Is has been shown that states which were 
obtained are true ground states or are at least very close 
to them. Even if one calculates only almost true ground 
states the resulting value of the stiffness energy is very 
reliable. It has been found that the stiffness energy in- 
creases with system size. So strong evidence has been 
obtained that the 3d Ising spin glass has a non-zero tran- 
sition temperature Tc. 

The only uncertainty arises from the fact that these 
calculations were restricted to systems sizes L < 10, it 
means anyway that the number of spins is more than 200 
times larger than for the systems examined before {L = 
6) 10-^. Additionally in contrast to former publications 
it is possible to estimate the quality of the low energy 
states. 
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TABLE L Simulation parameters: L — system size. Mi = 
initial size of population, v = average number of offsprings 
per configuration, rimin = number of CEA minimization steps 
per offspring, r = average computer time per ground state 
on a 80MHz PPC601, Nl ~ number of realizations of the 
random variables. 
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FIG. 1. Stiffness energy A as function of system size L. The 
line represents the function A{L) — aL^^ with Qs = 0.19(2). 
The inset shows the same figure on log-log scale. The increase 
of A with system size indicates, that in 3d Ising spin glasses 
an ordered phase exists below a non-zero temperature T^. 
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FIG. 3. Average stiffness energy A as function of simula- 
tion parameters Mi * u * rimin for 90 systems of size L = 6. 
This quantity is more sensitive to changes in the simulation 
parameters than the energy. 
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FIG. 2. Average energy as function of simulation parame- 
ters Mi *u * nmin for 90 systems of size L = &. With increas- 
ing numerical effort the energy lowers. For the two rightmost 
points [Mi — 64, 128; v — 4; rimin = 2) no further decrease of 
the energy is possible. A comparison with exact ground state 
calculations using a Branch-and-Cut algorithm confirms that 
in fact for all 90 systems the true ground states are found! 



FIG. 4. Average stiffness energy A as function of sample 
size Nl for three different sets of simulation parameters for 
system size L = 6. If the calculated states are close to the 
true ground state (set (16,4,2)) the resulting stiffness energy 
is very close to the true value and converges towards it with 
increasing sample size. (The 90 samples from Fig. ^ are not 
included in the samples here, so the values for Nl = 100 are 
different from the results above.) 
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